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Abstract 

Using numerical and analytic methods, we study the behavior of 
granular particles contained in a vibrating box. We measure, by molec- 
ular dynamics (MD) simulation, several quantities which characterize 
the system. These quantities — the density and the granular temper- 
ature fields, and the vertical expansion — obey scaling in the variable 
x = Af. Here, A and / are the amplitude and the frequency of the vi- 
bration. The behavior of these quantities is qualitatively different for 
small and large values of x. We also study the system using Navier- 
Stokes type equations developed by Haff. We develop a boundary con- 
dition for moving boundaries, and solve for the density and the tem- 
perature fields of the steady state in the quasi-incompressible limit, 
where the average separation between the particles is much smaller 
than the average diameter of the particles. The fields obtained from 
Haff 's equations show the same scaling as those from the simulations. 
The origin of the scaling can be easily understood. The behavior of 
the fields from the theory is consistent with the simulation data for 
small x, but they deviate significantly for large x. We argue that 
the deviation is due to the breakdown of the quasi-incompressibility 
condition for large x. 

PACS Number: 05.20.Dd, 46.10+z, 46.30.My, 61.20.Ja 
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1 Introduction 



Systems of granular particles (e.g. sand) exhibit many interesting phenom- 
ena, such as segregation under vibration or shear, density waves in the out- 
flow through a hopper and a tube, and the formation of heaps and convection 
cells under vibration |], ^ ||, These phenomena are consequences of the 
unusual dynamical response of the systems, many of which are still poorly 
understood. 

Granular particles in motion lose energy due to inelastic collisions of the 
particles. Thus energy has to be supplied to granular systems from an ex- 
ternal source (s) in order to sustain the movement of the particles. We can 
classify granular systems according to the way in which energy is supplied, or 
how the systems are agitated. A few common ways to agitate the systems are 
shear, vibration and body force (e.g. gravity). Here, we consider only vibra- 
tional agitation, which is probably one of the most popular method. There 
are many interesting phenomena associated with the vibrated systems, such 
as convection cells || ||, [7|, || , heap formation || [L(], [II], [12], |13| , subharmonic 



instability |H| , surface waves [TTS], |T(| and even turbulent flows |T7]]. The basis 
for understanding these diverse phenomena is, in our opinion, to understand 
of the state of granular media under vibration. The state is characterized by 
several fields of the system, for example, the density, velocity and granular 
temperature fields. 

There have been several studies on the state of granular systems under 
vibration. Thomas et al studied the system in three dimensions, mainly 
focusing on the behavior of shallow beds | I3fl . Clement and Rajchenbach 
experimentally measured the density, velocity and temperature fields of a 
two-dimensional vertical packing of beads |l9fl . They found that the tem- 
perature increases monotonically as the distance from the bottom plate is 
increased. The same system was studied by molecular dynamics (MD) sim- 
ulation with similar results 
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In a series of simulations and experiments, 
Luding et al studied the behavior of the one and two dimensional systems 
21 , |2~2 , p3| . They found that the vertical expansion, which is the increase of 
the center of mass due to the vibration, scales in the variable x = Af. Here, 
A and / are the amplitude and the frequency of the vibration. The expan- 
sion behaves as x 2 and x 3//2 in one and two dimensions, respectively. The 
reported shape of the temperature fields is decreasing as the distance from 
the bottom is increased, in contrast to the previous result |21|, [23| . In a recent 
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MD simulation of the three dimensional system, Lan and Rosato measured 



the density and temperature fields [24]. They compared the results with the 



theoretical predictions by Richman and Martin P5| |, and found good agree- 
ments. More discussion on the theory will come later. They also found the 
shape of the temperature field is different for low and high amplitude of the 
vibration. Finally, an approximate theory was developed for the system in 
one dimension, which agrees with simulations in the weak and the strong 



dissipative regime [26 



In this paper, we measure the density, the granular temperature and the 
expansion for the two dimensional system using MD simulation. We find that 
not only the expansion, but also the density and the temperature fields, scale 
in the variable x. The behavior of these quantities is qualitatively different 
for small and large values of x. The different shapes of temperature field 
reported can be due to different values of x used in the measurements. We 
also study the system using continuum equations developed by Haff |27[]. We 
develop a boundary condition for moving boundaries. We solve the equations 
for the density and the temperature fields and the expansion of the steady 
state. We consider only the quasi-incompressible case, where the average 
separation between the particles is much smaller than the average diameter of 
the particles. The fields obtained from Haff 's equations show the same scaling 
as those from the simulations. Furthermore, the origin of the scaling can be 
traced back to the boundary conditions, and can be easily understood. Also, 
behaviors of the fields from the theory are consistent with the simulational 
data for small x, but they deviate significantly for large x. We argue that the 
discrepancy is due to the breakdown of the quasi-incompressibility condition 
for large x. Some of these results can also be obtained using the theory 
by Richman and Martin, which is based on an extension of the Boltzman 



equation for inelastic gas [28, 29 



The paper is organized as follows. In Sec. [2], we start from defining the 
interaction of the particles used in the MD simulations. The geometry and 
various fields of the system will be discussed. We then present the expansion 
and the fields obtained by the simulations. Analytic results will be discussed 
in Sec. [| The continuum equations for granular material will be given, 
and simplified for the present system. We then derive a boundary condition 
for moving walls, and present the solution of the equations. In Sec. £|, we 
compare the fields obtained by the simulations with those by the theory. 
Finally, conclusions are given in Sec. R. 
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2 Molecular Dynamics Simulations 



2.1 Interaction of the Particles 

We start by describing the interactions used in the molecular dynamics sim- 
ulations. The force between two particles i and j, in contact with each other, 
is the following. Let the coordinates of the center of particle i (j) be Ri (Rj), 
and r = Ri — Rj. In two dimensions, we use a new coordinate system defined 
by two vectors h (normal) and s (shear). Here, n = f/\f\, and s is defined as 
rotating n clockwise by tt/2. The normal component i^jL»j of the force acting 
on particle i by j is 

Ff^ = k n {ai + aj - \r\) - j n m e (v-n), (1) 

where a« (aj) is the radius of particle % (j), and v = df/dt. The first term 
is an linear elastic force, where k n is the elastic constant of the material. 
The constant 7 n of the second term is the friction coefficient of a velocity 
dependent damping term, m e is the effective mass, mjTOj/(mj + rrij). The 
shear component F?_^ is given by 

F-^i = -%m e (v- s) - sign(5s) jmn{k a \5s\,fi\FJ_^\), (2) 

where the first term is a velocity dependent damping term similar to that 
of Eq. dH). The second term is to simulate static friction, which requires 



a finite amount of force (fiFji^) to break a contact [30]. Here, fj, is the 
friction coefficient, 5s the total shear displacement during a contact, and k s 
the elastic constant of a virtual spring. There are several studies on granular 
systems using the above type of interactions |31| . However, only a few of 



them [^, 31, |32j include static friction. A particle can also interact with a 
wall. The force on particle i, in contact with a wall, is given by Eqs. (0)-(@) 
with aj = oo and m e = A wall is assumed to be rigid, i.e. it is not moved 
by collisions with the particles. Also, the system is under a gravitational field 
g. We do not include the rotation of the particles in present simulation. A 



detailed explanation of the interaction is given elsewhere |33 . 

The trajectories of the particles are calculated using a fifth order predictor- 
corrector method. The interaction parameters used in this study are fixed as 
follows, unless otherwise specified. They are k n = 5 x 10 6 , k s = 1 x 10 4 ,7„ = 
1 x 10 3 , 7 S = and \l = 0.2. The timestep is taken to be 5 x 10~ 6 . This small 
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timestep is necessary for the large elastic constant used in the simulations. 
For too small values of the elastic constant, the system loses the character of 
a system of distinct particles, and behaves like a viscous material. In order to 
avoid artifacts of monodisperse systems (e.g., hexagonal packing), we choose 
the radius of the particles from a Gaussian distribution with the mean 0.1 
and the width 0.02. The density of the particles is 0.5. Throughout this 
paper, CGS units are implied. 



2.2 Setup for the Simulations 

We put the particles in a two-dimensional rectangular box. The box consists 
of two horizontal (top and bottom) plates which oscillate sinusoidally along 
the vertical direction with a given frequency and amplitude. The separation 
between the two plates is chosen to be much larger (10 5 times) than the 
average radius of the particles, so the particles never interact with the top 
plate for all the cases studied here. For simplicity, we apply a periodic 
boundary condition in the horizontal direction. 

We start simulation by inserting the particles at random positions in the 
box. We let them fall by gravity and wait until they lose energy by collisions. 
We wait for 10 5 iterations for the particles to relax, and during this period we 
keep the plates fixed. The typical velocity at the end of the relaxation is of 
order of 10~ 2 . After the relaxation, we start oscillating the plates. We vibrate 
for about 50 cycles before taking any measurements in order to eliminate any 
transient effect. 

The main quantities we are interested in are density, velocity and granular 
temperature, which probably are the most basic quantities to characterize 
the system. In order to measure these fields, we divide the system into 
rectangular cells of width w and height h. For cell i, we identify the particles 
whose center lies within the cell, which we label by j = 1,2,3..., iVj. Here, 
Ni is the total number of particles in the cell. The density at cell i, p iy is 
defined as 

I ^ 

Pi = r (3) 

w ■ h p{ J 

where rj is the radius of particle j. We define an average of a certain quantity 
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A in cell i, (A)i, as 

(A), = (4) 

where m,j is the mass of particle j. Then, the velocity field at cell i,Vi, is 
given by 

Vi = (v)i, (5) 
and the granular temperature at cell i, T { is 

Ti=\{{vl-{v x )% + {vl-{v y )%). (6) 

We also measure the vertical component of the center of mass y cm of the 
particles, which is an useful quantity to compare with theoretical predictions. 
In principle, we can calculate y cm from the measured density fields, but the 
calculation introduces additional sources of uncertainty. We measure the 
above quantities by averaging over at least three independent runs, where 
each run is averaged over about 200 cycles of vibrations. 

We study the effects of the width W of the vibrating plates on the mea- 
sured quantities in order to determine the optimal width for the main mea- 
surements. We measure the density and the granular temperature fields for 
W = 1,2,3,4, where the height H of the pile is roughly fixed constant by 
changing the number of particles to be 50, 100, 150, 200. The resulting height 
is little smaller than 2. Also, we use w = W and h = 0.2. Here, the large 
w is chosen due to the translational symmetry resulted from the periodic 
boundary condition. We do not find any systematic dependence of the fields 
on W except that they become less noisy for larger W, which is simply due 
to the larger number of particles contained in a cell. We find W = 1 is a 
good compromise between the quality of data and the computation time. 
The dependence of the fields on H will be discussed later. 



2.3 Measurement of the Expansion 

We now discuss the quantities obtained as described previously. We start 
from the vertical expansion of the pile y exp , defined as the difference in the 
vertical center of mass y cm during and before vibration. Since the expansion 
can be expressed as an integral involving the density field, y exp can be used as 
a number representative of the density field. In Fig. 1, we show y exp measured 
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for various values of the amplitudes A and the frequencies / of the vibration 



34}|. Here, the width W = 1 and the number of particles N = 50. Note that 
the expansion is plotted against a scaling variable x = Af. For the entire 
range of variable x, we do not find any systematic deviation from the scaling 
behavior except that the lowest frequency data (/ = 20) seems to be slightly 
off. This scaling behavior was first measured in the simulations by Luding 



et al [Ell 22, 23]. For intermediate values of x, our data is also consistent 



with the x 3 / 2 behavior proposed by Luding et al |22] . However, for small and 



large values of x, our data shows deviations from the behavior. The pattern 
of the deviation is not affected by changing the values of the elastic constant 
k n (to 10 6 and 5 x 10 5 ). The scaling behavior and the deviation will later be 
discussed in more detail. 

Since we are concerned about the possibility that y cm is dominated by a 
few top particles isolated from the main part of the pile, we also try an alter- 
native definition of the expansion. We sort the particles in a box according 
to their vertical coordinates. We then define y ce nter as the vertical coordinate 
of the center of the particle which is the (iV/2)-th in the list. Here, N is the 
total number of particles in the box. The quantity y ce nter still contains the 
meaning of the center of mass, but is now not dominated by a few particles. 
We find that, in contrast to our worries, the expansions using y center behave 
essentially the same way as those using y cm . We thus use the definition using 
y cm from now on. 



2.4 Measurement of Density and Temperature 

We discuss the density p(y) and the granular temperature T(y) fields. First, 
we study the behavior of the fields at the bottom plate, i.e., in the lowest 
cell (y — 0). These values are, as we shall see later, easier to compare with 
theoretical predictions than the entire fields. In Fig. 2 (a), we show the 
density p(y = 0) vs x for various values of A and /. Here, we find the 
same scaling behavior as the expansion y exp (Fig. 1). This is not surprising 
considering the relation between the expansion and the density field. We 
also find the same scaling behavior for the temperature. The square root 
of T(y = 0) is plotted against the scaling variable x in Fig. 2(b). Although 
there are larger fluctuations, the same scaling is still apparent. Also, the 
decay of T(y = 0) for large x is partly due to a dip in the density near the 
bottom plate. 
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We then proceed to study the entire ^/-dependent fields. In Fig. 3, we 
show the density field p(y) for several values of A and /. Here, the value of 
x is fixed to be 3 and 10 in Fig. 3(a) and (b), respectively. In the figures, 
the quality of scaling is poorer than those shown previously (Figs. 1 and 2). 
However, the quality of scaling will be improved by discarding the lowest 
frequency data (/ = 20). The density field plotted against y/y exp was also 
shown to collapse into a single curve p3| , which is consistent with above 



scaling, since y exp scales in x. The situation is not much different for the 
temperature field T(y). In Fig. 4(a) and (b), we show the square root of 
T(y), where x is again fixed to be 3 and 10, respectively. Here also, the 
quality of scaling is not great, and will be improved by ignoring the / = 20 
data. Putting aside the deviations at low frequencies, it is not entirely clear 
from these data that either the density or the temperature field exhibits more 
than an approximate scaling. 

One interesting point is the shape of the temperature field. For large 
values of x, the field T(y) is a monotonically decaying function of the height 
y (as shown in Fig. 4(b)). As one decreases x, a local maximum begins to 
appear. In Fig. 4(a) one can see an maximum is starting to form around 
y = 2. As one decreases x further, the maximum becomes more and more 
pronounced. The above observation can explain the different shapes of T(y) 
reported by several experiments and simulations on the granular systems. 
In experiments of the two dimensional system, a local maximum is always 



present in the T(y) field []19 |. In contrast, T(y) is found to be monotonically 
decaying in simulations of the one dimensional systems [21, 2l|. Other simu- 
lations of the system in three dimensions also report monotonic decay of the 
temperature field. In the three dimensional simulations, however, changes of 
the shape of the field for smaller amplitudes were noted, although were not 
systematically studied [p4"| . 

We conclude this section by presenting the dependence of y exp on the total 
height H of the pile. We change the total height by varying the number of 
particles N and keeping the width W fixed. We keep the previous value of 
W = 1. In Fig. 5, we show the scaled expansion for N = 50, 100, 150 and 
for several values of A, where / is fixed to 100. The expansion is scaled by 
50/ N, where the factor 50 is used to keep the N = 50 data unchanged. In 
the figure, we find a reasonable collapse of the curves. There seems to be no 
systematic deviation, although the data for large x seems to be slightly off. 



The scaling is first found by Luding et al [21 
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In this section, we have presented the data from the MD simulations. 
The expansion as well as the density and the granular temperature fields 
exhibit scaling in the variable x, but the quality of the scaling for the fields is 
worse than that of the expansion. Also, the lowest frequency data seems to 
deviate from the scaling. In the next section, we present a continuum theory 
of granular material which can explain the origin of the observed scaling. 



3 Analytic Results 
3.1 Basic Equations 

The system of granular particles needs an external agitation(s) in order to 
sustain the movement of the particles. Without such an agitation, the motion 
of the particles will be suppressed because of the constant energy loss result- 
ing from inelastic collisions of the particles. Many granular systems with an 
external agitation, including the system of granular particles in a vibrating 
box, can loosely be described as inelastic gas. A particle in such a system con- 
stantly moves. It assumes a ballistic trajectory until it collides with another 
particle or, if there is, a wall of the container. After the collision, the particle 
again moves ballistically. Such movement is analogous to those of particles 
in a gas. Motivated by the analogy, tools used to develop the kinetic theory 
of gas have been adapted to study granular systems. One such effort is to 
extend the Boltzman equation to an inelastic gas |58, |29|. In this approach, 



the time evolution of the fields (velocity, density and granular temperature) 
is written in terms of an integral involving all possible collisions, namely the 
Boltzman integral. Another approach, due to Haff, is to develop equations 
of motion very similar to the Navier-Stokes equations |27|. However, the 



coefficients in the equations, e.g. viscosity, are no longer constant, but now 
functions of the fields. Both approaches are not yet complete, and they con- 
tain assumptions whose validities have yet to be checked. Haff's approach 
seems to be less rigorous; for example, the dependence of the coefficients on 
the fields are derived based only on intuitive arguments. On the other hand, 
the equations one has to solve in Haff's approach are much simpler. Despite 
its simplicity, known solutions for a few systems using Haff's approach do 
not significantly differ from those using the other method |27|, ^9|. In this 
paper, we will follow Haff's approach. 
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The equations of motion developed by Haff consist of mass, momentum 
and energy conservation. The mass conservation equation is 

lp + V • (pv) = 0, (7) 

where p and v are the density and the velocity fields, respectively. The 
equation (|7p is exactly the same as that of the Navier-Stokes equations. Next 
is the i-ih component of the momentum conservation equation, 

d -> d -> Ov ■ Ov ■ 

P - Vi + p( ff. V)v, = ^[-p + A(V • 8)] + + + ». (8) 

where summation over index j is implied. The coefficients A and 77 are vis- 
cosities which will be determined later. Also, p is the internal pressure, and 
is the i-th component of the gravitational field. The form of Eq. (§) is again 
the same as that of the Navier-Stokes equations. However, the coefficients as 
well as the internal pressure are now functions of the fields instead of being 
constant. The last of the equations of motion is energy conservation, 

^-(-pv 2 + -pT) + J^[(-pv 2 + -pT)v i ] (9) 



dx. 
+ pVi9i 

Here, T is the granular temperature field defined in Eq. @, K is the "ther- 
mal conductivity," and I is the rate of dissipation due to inelastic collisions 
35f] . Also, the summation over repeated indices is implied. Although the 



form of Eq. (|^) is somewhat different from that of the Navier-Stokes equa- 
tions, the equation can still be easily understood. The left hand side of Eq. 
@j is simply the material derivative of the total kinetic energy, where the 
total kinetic energy is divided into the convective part (involving v) and the 
fluctuating part (involving T). On the right hand side of the equation, the 
first three lines are simply the rate of work done by the internal pressure, 
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viscosity and gravity, respectively The term involving K is the rate of energy 
transported by "thermal conduction." The term /, which is a consequence 
of inelasticity of the particles, is responsible for many unique properties of 
granular material. 

We now discuss the above coefficients which are yet to be determined. 
The relations of these coefficients to the fields are derived based on intuitive 
arguments. Also, the derivation assumes that the density is not significantly 
smaller than the closed packed density, i.e., the system is almost incompress- 
ible. The relation for the internal pressure is 

T 

p = tdp—, (10) 

s 

where t is an undetermined constant, d is the average diameter of the par- 
ticles. The variable s, which is roughly the gap between the particles, is 
related to the density by 

where m is the average mass of the particles. Then, the viscosity 77 is given 

as 

v = qd 2 P —, (12) 

where q is an undetermined constant. In a similar way, the thermal conduc- 
tivity is found to be 

K = rd 2 — . (13) 

s 

Here again, r is an undetermined constant. Finally, the rate of dissipation is 

rp3/2 

I = 19 • (14) 

s 

Also, 7 is an undetermined constant. One can notice the viscosity A is 
left undetermined. This is due to the fact that in the range where these 
relations are valid, the term containing A is negligible, and is dropped from 
the calculation. 
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3.2 Equations for Particles in a Vibrating Box 

We now apply Haff 's equations to the system being considered in the present 
paper — granular particles under vibration. Since the general solution of the 
problem using Haff's equations are too difficult to obtain analytically, we 
introduce several constraints which simplify the equations. First, we only 
study steady state properties of the system. In order words, we are interested 
in only the time averaged quantities. Since there is no net flow of particles in 
the steady state, the time averaged values of v is zero, which greatly simplifies 
the equations. The second simplification is resulting from the horizontal 
periodic boundary conditions we imposed. Due to the boundary condition, 
there is no significant variation of the fields along the horizontal direction. 
As a result we only have to deal with a one dimensional equation instead 
of two or three dimensional ones. The third condition is incompressibility, 
which is a little bit tricky. Incompressibility implies, strictly speaking, that 
the density p is constant. Due to the relation between p and s (Eq. (|TT|)), s 
also has to be constant. Here, we are interested in the situation where s is 
much smaller than d, but still non-zero. In such cases, the variation of the 
density can be ignored, but not the variations of the variables that depends 
directly on s. We call this condition quasi-incompressibility. 

The simplified equations for the pile of granular particles under vibration 
are then 

-, = 0, (15) 
-^p + pg = 0, (16) 

and 

I p A[^i_ T ]_/ = o. (17) 
Tdy V dy 1 K ' 

Here, all the fields are averaged over time, so they are not functions of time, 
but of spatial coordinate y only. We thus replace partial derivatives in the 
equations with ordinary derivatives. Note that the first equation is automat- 
ically satisfied. Substituting the relation for the internal pressure Eq. fllTf ) to 
Eq. (jlB]), we obtain 

£<7> = -& (18) 
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whose solution is 

7 = ->-*). as) 

where y Q is a constant which will be discussed later. Substituting Eq. (|T^) 
to Eq. ( pj) , we obtain 

-f^VT + -^VT - ^-VT = (20) 

ay 2 V — Vo dy r d 



which is a linear differential equation for yT. We make a change of variable, 

z= V -^, (21) 



where t is defined as cf. By inserting the new variable to Eq. (|20|) , we 

arrive at the final form of the equation 

d^ t — 1 d i — i — . . 

— VT+-— Vf- Vf = 0. 22 

az A z az 

The solutions the equation are the modified Bessel functions I {z) and K Q (z). 
However, K Q (z) has a singularity at the physical region, and can not be a 
part of the solution. The solution is, therefore, 

VT=AI (z), (23) 

where A is a constant yet to be determined. The result obtained so far is 
identical to that of Haff |27| . We now discuss our contribution to the problem, 
which concerns the boundary conditions. 

3.3 Boundary Condition 

One of the areas which are yet to be completed in Haff 's theory is boundary 
conditions. For example, what is the boundary condition one has to impose 
for the temperature field T and the velocity field v at the bottom of a vi- 
brating box? The boundary conditions for fixed walls have been developed 
based on energy conservation across the boundary |36|. We will show that 



a boundary condition for moving walls can also be derived from the similar 



condition of energy conservation. As noted in [36], there are two energy flows 
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near a boundary. One is the thermal conduction, and the other is the energy 
transfer by collisions between the particles and the boundary. Energy con- 
servation implies that the rate of energy transferred by the conduction has 
to be equal to that by the collisions. We consider an horizontal boundary 
vibrating vertically, like the bottom plate of a vibrating box. The rate of the 
thermal conduction per area at the boundary is given as 

\ K % T I- < 24 > 

where y = is the average vertical coordinate of the boundary. To calculate 
the rate of energy transferred by the collisions, consider a particle colliding 
with the boundary. If convection is absent, the movement of the particle 
is resulted from the fluctuation. The typical velocity due to the fluctuation 
is y/T. The velocity of the particle after the collision will depend on the 
velocity of the boundary at the moment of the collision. Here, we assume 
that the phase of the vibration at which the collision occurs is random — an 
assumption whose validity will be discussed later. In such cases, the typical 
velocity of the wall is v w = Af. The amount of energy transferred in one 
collision is 

\m[T(l-el)-(l + e w fvll (25) 

where e w is the coefficient of restitution between a particle and the wall. 
The rate of energy transfer per area is Eq. (|25|) multiplied by the frequency 
of collisions T/s divided by the area of contact d 2 . Here, we again use the 
quasi- incompressibility. Equating this rate to Eq. ([24]) , and use the definition 
of K (Eq. (|I3|)), we obtain 

T{l-el)-{l + e w fvl = -^ (26) 

a ay 

where a is an undetermined constant. The method used to derive the above 
boundary condition can be extended to a case where the boundary is moving 
horizontally, and even to a case with both horizontal and vertical motions. 

We now apply the above boundary condition to the problem of the vi- 
brating box. First, note that the modified Bessel function I {z) can be ap- 
proximated as 

«*) ~ (27) 
\/2irz 
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for z > 0. The expression is actually an asymptotic form of I {z) for z ^> 1, 
but it still is a good approximation even for small values of z. By applying the 
boundary condition Eq. fl2"6| ) to the solution Eq. (f23|) with the approximation 
Eq. fl27|), we obtain 

/ 2 (-W^) l-el-(2rd/a£) (l-£/2y ) w 1 J 



' B 2 v 2 



w '■ 



where a new variable B is introduced. Substituting for A in the solution 
Eq. fl23|), we now arrive at the final form of the solution for the temperature 
field , 

y 



T(y) ~Bv w J-?— exp(-y/e), (29) 

\Va-v 

where we use the property I (— z) = I (z). Besides the prefactor B, the 
above result is the same as that obtained by Haff who simply assumed that 



'T at the boundary is v w |27|]. The temperature field decays exponentially 
with decay length I as one moves away from the bottom plate. Note that 
it has a singularity at y = y a , which is an artifact resulting from the fact 
that Haff 's theory is not valid near the top part of the pile. In other words, 
the theory is based on the picture that the state of the system is gaseous. 
Near the top of the pile, however, the particles do not collide frequently, but 



assume ballistic trajectories p7[ . On the other hand, Eq. ( p23|) is still valid 
for y not very close to y Q , since the gaseous picture is valid for that ranges of 
y. Thus, the behavior of T to increase again for large y is a physical effect, 
not an artifact. A word about y a . One can see y Q is still left undetermined. 
According to Eq. (0), the condition to determine y a is that the internal 
pressure vanishes at y — y . However, since the solution we have obtained is 
not valid at y = y , we can not impose the condition to determine y . 

Having determined T(y), we now determine the density field, represented 
here by the separation s(y). Inserting the solution of T(y) into Eq. (|T9"|), we 
obtain 

s(y) * ^B'vl^-^eM-^y/i), (30) 

where p Q is the density of the pile at rest. The separation s also decays 
exponentially as one moves away from the bottom plate, then increases for 
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large values of y, similar to T(y). Also, the above expression for s is not valid 
in the neighborhood of y — y Q . Finally, we consider the expansion y exp . The 
expansion can be written as 



4 Discussion 

4.1 Condition for the Gaseous State 

In the previous section, we have obtained analytic results for the vibrating 
box problem using Haff 's theory. There are several assumptions made in the 
construction of the theory as well as in the derivation of the solution. In this 
section, we check the validity of the solution, and in turn the assumptions, 
by comparing with the results obtained by the MD simulations. 

Before making the comparison, we check the condition under which the 
system is regarded as gaseous. We consider that the system is gaseous if the 
particles in the system stay afloat for the most part of time colliding with 
each other. Therefore, for the gaseous system, the vibration should provide 
enough input to make the particles float. In other words, the collisional force 
caused by the bottom plate should be large enough to overcome the pressure 
due to the gravity. Consider a particle at the bottom of the box, where the 
pressure is the greatest. The typical change of momentum during a collision 
of the particle with the bottom plate is mv m , where m is the mass of the 
particle, and v w is the typical velocity of the plate. Here, we choose v w to be 
27rAf, the maximum velocity of the plate. We then consider the frequency 
of such collisions. We argue that the particle collides roughly once during 
one period, at least for the small amplitudes of the vibration. It will be later 
shown that the system becomes gaseous at a small value of the amplitude, 
so the above estimation seems to be reasonable. The total force generated 
due to the collisions is of order of 27rmAf 2 . On the other hand, the force on 
the particle due to gravity is of the order of mgn, where n is the number of 
layers in the pile. Then, the condition for the gaseous system becomes that 
the ratio F of the collisional force to the gravitational force to be much larger 




(31) 



t 2 dp 
2gh Q 
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than one. The ratio F is essentially the Froude number. The condition can 
be written as 

F ~ » 1, (32) 
Ann 

where T is 4ir 2 Af 2 /g, the maximum acceleration of the plate scaled by that 
of the gravity. Thus, the acceleration of the bottom plate, not the velocity, 
is relevant in determining the condition for the system to be gaseous. 

We check the condition also by MD simulation. At a given time, we 
measure the total number of pairs of particles in contact with each other. We 
average the number over 200 periods of the vibration. The system is gaseous, 
if the averaged number of pairs N p is much smaller than the total number of 
particles N. Here, we use the criterion that N p < 0.1 AT for the system to be 
gaseous. We choose the width W = 1 and N = 50. For several frequencies, 
we increase the amplitude A until N p /N becomes less than 0.1. We obtain, 
for each /, an interval of A containing the value at which N p /N = 0.1. 
The intervals are [0.06, 0.08], [0.01, 0.02] and [0.003, 0.004] for / = 20, 50, 100, 
respectively. The corresponding intervals of V are [0.97, 1.29], [1.01, 2.01] and 
[1.21, 1.61]. The acceleration at which the system becomes gaseous is roughly 
constant, which is consistent with the above argument. We now study the 
height dependence of the condition. We fix the frequency / = 100 and the 
width W = 1, and vary iV to 100 and 150. By using the same criterion, 
we find the intervals of T to be [2.41,3.22] and [4.03,5.64] for N = 100 
and 150, respectively. The acceleration needed for the gaseous system is 
consistent with a linear increase in N, which is again in agreement with the 
above argument. One thing to note is that T to make the system gaseous is 
fairly close 1 for iV = 50. For the values of V used in the simulations of the 
preceding section, the system is always gaseous. Also, N p for large amplitudes 
A (~ 0.1) is essentially zero, which will be further discussed later. 



4.2 Scaling of Density and Temperature 

We compare the expansion y exp measured by the MD simulations with that 
predicted by the theory. Measured values of y exp , as shown in Fig. 1, exhibit 
good scaling for the variable x = Af with no systematic deviation for the 
entire range of x. The scaling is exactly what is predicted by the theory 
in Eq. fl3~ll) . The theory also predicts y exp is inversely proportional to the 
rest height of the pile H for fixed A and /, which is certainly consistent 
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with the data presented in Fig. 5. On the other hand, y exp should behave 
as x 2 according to the theory, but the data (Fig. 1) is consistent with the 
behavior only for small values of x (up to ~ 2). For larger values of x, y exp 
definitely increases slower than x 2 . The situations for the density p(y) and 
granular temperature T(y) fields are very similar. The value of the fields at 
the bottom plate (y = 0) shows good scaling for x (Fig. 2), as again predicted 
in Eqs. (29) and (|30|) . The value of T(0) from the simulations increases for 



small x, reaches a maximum around 4, than decays, very different from the 
x 2 behavior given by Eq. (P9|) . The behavior is again consistent with the 
prediction only for small values of x. For the density field, we recall the 
relation between the separation s and the density p, which is 

where p Q is the density of the rest pile. In Fig. 6, we plot p(0)~ 1//3 against 
x 2 . The resulting curve, as predicted in Eq. (j30"l) , should be straight. The 
curve is indeed straight, but only for small values of x. For large x, the curve 
increases slower than linearly. 

All three quantities we have studied behave in a similar manner. They 
all scale in the variable x, they are all consistent with the theory for small 
values of x, and they all deviates from the predicted behavior as x is increased. 
Furthermore, the value of x where they start to deviate seems to be roughly 
the same (x = 2 ~ 3). We recall that one of the assumptions in deriving the 
solution in the previous section is quasi-incompressibility, which implies that 
the separation s is much smaller than d, the average diameter of the particles. 
However, as shown in Fig. 2(a), the density of the pile becomes the half of a 
packed density for x ~ 3, which is in violation of the assumption. The system 
can not be treated as quasi-incompressible for values of x larger than about 
3. Since the point where the quasi-incompressbility breaks down seems to 
coincides to the point where the deviation from the theory starts, we strongly 
suspect that the deviations are resulted from the breakdown. Furthermore, 
if the deviation is due to the breakdown, it is clear why the relevant variable 
for the breakdown is x, not T or other variable. The breakdown occurs when 
s is comparable to d, and s scales in the the variable x. Thus, the breakdown 
occurs at a specific value of x. 

We now consider the behavior of the entire density and temperature fields. 
The density field p(y), as shown in Fig. 3, roughly scales in x and has one 
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maximum, both of which are consistent with Eq. (|30|). The temperature 
field T(y) also scales in x, but displays qualitatively different shapes (Fig. 4) 
depending on A. For small A, the field has a local maximum around y = 

2, while the maximum disappears for larger A. On the other hand, the 
maximum is always present in the prediction of the theory Eq. (|29|) . We study 
the condition for the change of the shape. For given /, we increase A until the 
local maximum disappears. We obtain intervals of A containing the value 
at which the maximum disappears. They are [0.10, 0.15], [0.04, 0.06] and 
[0.019, 0.020] for / = 20, 50, 100, respectively. The corresponding intervals of 
x for the different frequencies are roughly the same. We thus face a familiar 
situation — the behavior of T(y) is consistent with the theory for small x, but 
deviates for large x. Furthermore, the value of x at which T(y) deviates from 
the theory is about 2, which is again consistent with the previous argument. 

4.3 Beyond Incompressibility 

Since we argue quasi-incompressibility is valid only for modest range of x, 
one might wonder why we do not relax the condition. The simple answer 
is that the system of equations becomes too complicated. We now consider 
the changes needed in order to relax the quasi-incompressibility. The three 
equations Eqs. (H)-® are written in general form, and do not need any 
modification. The relations of p,r/,K,I to the fields Eqs. (|10|)-(1I4|) as well 
as the boundary conditions Eq. (^) have to be modified. It is not these 
modifications themselves, but the complexity of the resulting equations, that 
makes the analytic solution too difficult to obtain. A word about the regime 
of validity. The range of x where the quasi-incompressible theory is valid 
depends on H the rest height of the pile. For H = 2, which corresponds 
to the system of about 10 layers, the theory is valid until x reaches about 

3. As we increase H, the range of validity also increases. Since a typical 
experiment of granular system involves several tens or more layers, the quasi- 
incompressible theory can describe quite large range of x, and we hope many 
interesting phenomena occur within the regime of validity. 

One also might wonder why the scaling predicted by the theory is valid 
even in the regime where the quasi- incompressible theory is no longer valid. 
The only place where the external parameters A and / come to the system 
is the boundary condition. The boundary condition, as previously discussed, 
is a consequence of energy conservation. Even in the general situation of 
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a compressible system, the boundary condition is obtained by equating two 
energy fluxes. One is the energy transferred by collisions between the parti- 
cles and the wall, which is a function of the velocity of the wall. Now, the 
velocity is proportional to Af, but the other flux term does not depend on 
A or /, but on the fields of the system. Thus, it is not surprising to see the 
scaling for the variable Af holds even for the compressible systems. 

We finish the section by discussing two behaviors previously mentioned. 
As noted in Sec. |2]3|, the expansion y exp seems to increase faster near the 
largest value of x. This may be explained as follows. As A is increased for 
fixed /, the system becomes more expanded, and the number of interparticle 
collisions decreases. For sufficiently large A, it is be possible that a particle 
in the system hardly interact with other particles, and it mainly interacts 
with the bottom plate |]37|]. This is supported by the measurement of the 



average number of pairs N p in contact, discussed in Sec. pTT] . We find that 
N p is essentially zero near the largest values of A. If we completely ignore 
the interparticle interaction, it is straightforward to show that y exp ~ x 2 . 
Thus, y exp has to increase faster for large x, which is consistent with the data 
(Fig. 1). However, it is not clear, without further detailed study, that above 
mechanism is actually operating. Also, note that the mechanism predict 
scaling for variable x. The behavior of y exp was shown to depend on properties 
of the sidewall p2| , p3| . It was shown that y exp behaves as x 3 ^ 2 for entire range 



of x with inelastic sidewalls, while it deviate from the behavior with elastic 
ones (no energy loss). 

The other problem I want to discuss is the deviation from the scaling 
for very small values of /. In Sees. |273| and ^L4\ , we find nagging deviations 
from the scaling of the expansion, the density and the temperature field for 
/ = 20 data. One possible way to explain them is to consider the way the 
boundary condition is imposed. In constructing the boundary condition, we 
assume that particles collide with the bottom plate at random phase of the 
vibration. We thus assume the "thermal motion" is much larger than the 
convective motion. For low enough frequencies, however, the particles have 
enough time to dissipate their thermal fluctuations, so the assumption may 
not be valid. In such cases, the convective motion of the particle dominates, 
and one can no longer think of the system as particles with random motion. 
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5 Conclusion 



We study the behavior of granular particles contained in a vibrating box. 
Using MD simulation, we show that several quantities which describe the 
system — the density and the temperature fields as well as the expansion — 
obey scaling in the variable x = Af. The behavior of these quantities is 
qualitatively different for low and high values of x. We also study the system 
using continuum equations developed by Haff. We develop a boundary con- 
dition for moving boundaries, and apply at the bottom plate of the box. We 
solve for the density and the temperature fields of the time averaged steady 
state of the system. Here, we are limited to the quasi-incompressible case, 
where the average separation between the particles is much smaller than the 
average diameter of the particles. The fields obtained from Haff's equations 
show the same scaling as the simulational counterparts. The behaviors of 
the fields from the theory are consistent with the simulational data for small 
x, but they deviate significantly for large x. We argue that the discrepancy 
is due to the fact that the quasi-incompressibility condition we impose is not 
valid for large x. 

In this paper, we have showed that Haff's theory, even with many assump- 
tions, seems to describe the system reasonably well in the quasi-incompressible 
case (small x). It goes without saying that the current study is only a very 
small step towards the understanding of the system, and the following are 
what we consider important things yet to be understood. The quality of the 
simulational data needs to be improved in order to convincingly demonstrate 
the scaling. For example, it is not certain from Figs. 3 and 4 that p(y) and 
T(y) obey strict scaling or an approximate one. Also, the improvement is 
necessary for quantitative comparison of the fields from the theory with from 
simulations. Next is the validity of Haff's theory in the compressible regime 
(large x). As discussed in Sec. fO| , due to the complexity of the problem 
in the regime, it may be too difficult to obtain a closed form analytic solu- 
tion of the problem. However, with perturbative or approximate methods, 
one might be able to get some information. Also, one can numerically solve 
Haff's equations, and compare the results with MD simulations. Another 
issue to be understood are time dependent properties of the system. In this 
paper, we have studied only steady state properties, which averages out the 
variations of the fields during one period of the vibration. Time dependent 
properties can be studied by perturbation from the steady state p8| . And 
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of course, Haff 's theory has to be applied to other systems to find the power 
and limitation of the theory. 

The author thanks Joel Koplik and Stefan Luding for many useful dis- 
cussions and comments on the manuscript. This work is supported in part 
by the Department of Energy under grant DE-FG02-93-ER14327. 
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Figure Captions 



Fig. 1: The expansion y exp for several values of the amplitudes A and the 
frequencies / of the vibration. We find the expansion exhibits good 
scaling behavior with variable Af. 

Fig. 2: (a) The density and (b) the square root of granular temperature at 
the bottom plate for various values of A and / of the vibration. Here, 
we find that the same scaling as the expansion holds. 

Fig. 3: The density field for various values of A and /. Here, value of Af 
is fixed to be (a) 3 and (b) 10. We find that the same scaling as the 
expansion y exp seems to hold. However, the quality of scaling is poorer. 

Fig. 4: The square root of the granular temperature field for various values 
of A and /. The value of Af is fixed to be (a) 3 and (b) 10. We again 
find poorer scaling similar to the density field (Fig. 3). 

Fig. 5: The scaled expansion vs Af for N = 50, 100, 150 and for several 
values of the amplitude. All the curves seems to collapse to a single 
curve. 

Fig. 6: Plot of p(0) -1 / 3 against x 2 calculated from Fig. 2(a). The curve is 
straight for small values of x, and deviates from the linear behavior for 
large x. 
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